
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
    clear;
    clc;
    setEnv();
    workingDir = 'C:\Users\rniehaus\Documents\Elegans\';
    studyInstanceName = 'Study-RFUM294_loops';
    [dataList, suppData] = getStudyInstance(workingDir, studyInstanceName);

    for iFrame = 1:length(dataList)
        bwImage = suppData(iFrame).BwImage;
        
          %Plot bw image
             figure(1), imshow(bwImage) 
        
        %% GEt outer boundary
        [B,L,N,A] = bwboundaries(bwImage);

        outerBoundary = B{1};

        %Convert from subscripts to coordinates
        vidRows = size(bwImage,1);
        yOuter = vidRows - outerBoundary(:,1);
        xOuter = outerBoundary(:,2);
        outer = [xOuter,yOuter];
        
        %% Get inner boundary
            
        if length(B) > 1;
            innerBoundaryCell = B(2);
            innerBoundaryMat = innerBoundaryCell{1};
            yInner = vidRows - innerBoundaryMat(:,1);
            xInner = innerBoundaryMat(:,2);
            inner = [xInner, yInner];
             totalBoundary = [[xInner;xOuter],[yInner;yOuter] ];
        else
            xInner = [];
            yInner = [];
            inner = [];
            totalBoundary = [xOuter,yOuter ];
        end
        
              
            
        
            
        
        %% Get side 1 and side 2

        %Get Sides A (from head to tail) and Side B(from head to tail)of
        %the outer boundary
        
        
        if dataList(iFrame).Curv_EndARow == 1 && dataList(iFrame).Curv_EndACol == 1
            %Delta configuration
            endA = [dataList(iFrame).Curv_EndBCol, vidRows - dataList(iFrame).Curv_EndBRow ];
            dist = round(size(outer,1)/2);
            endB = [outer(dist,1), outer(dist,2)];
        elseif dataList(iFrame).Curv_EndBRow == 1 && dataList(iFrame).Curv_EndBCol == 1
            %Delta configuration
            endA = [dataList(iFrame).Curv_EndACol, vidRows - dataList(iFrame).Curv_EndARow ];
            dist = round(size(outer,1)/2);
            endB = [outer(dist,1), outer(dist,2)];
        else    
          
            endA = [dataList(iFrame).Curv_EndACol, vidRows - dataList(iFrame).Curv_EndARow ];
            endB = [dataList(iFrame).Curv_EndBCol, vidRows - dataList(iFrame).Curv_EndBRow ];
        end 
        
        %Find locations in boundary
        [~, locA] = ismember(endA, outer, 'rows');
        [~, locB] = ismember(endB, outer, 'rows');
        
        side1 = [];
        side2 = [];

        if locA < locB
            side1 = outer(locA:locB,:);
            side2 = [outer(locB:size(outer,1),:); outer(1:locA,:)];
        else
            side1 = outer(locB:locA,:);
            side2 = [outer(locA:size(outer,1),:); outer(1:locB,:)];
        end
       
            

        %Plot inner and outer boundary with sides shown
        figure(2)
        plot(side1(:,1), side1(:,2), 'g');
        hold on
        plot(side2(:,1), side2(:,2), 'm');
        
        if ~isempty(xInner)
            hold on
            plot(xInner, yInner, 'b'); 
        end
        hold off
            
         %% Delaunay

        tri = delaunay(totalBoundary(:,1), totalBoundary(:,2));

        %Plot the entire boundary and mesh
        figure(3)
        triplot(tri, totalBoundary(:,1), totalBoundary(:,2));    
            

            %% Create triangle structure        
             triangle = struct();
             for i = 1:size(tri,1)
                v1 = totalBoundary(tri(i,1),:);
                v2 = totalBoundary(tri(i,2),:);
                v3 = totalBoundary(tri(i,3),:);
                triangle(i).Vertex1 = v1;
                triangle(i).Vertex2 = v2;
                triangle(i).Vertex3 = v3;
                triangle(i).Vref1 = tri(i,1);
                triangle(i).Vref2 = tri(i,2);
                 triangle(i).Vref3 = tri(i,3);
                V12 = [[v2(1) - v1(1), v2(2) - v1(2)],0];
                V21 =  - V12;
                V13 = [[v3(1) - v1(1), v3(2) - v1(2)],0];
                V31= - V13;
                V23 = [[v3(1) - v2(1), v3(2) - v2(2)],0];
                V32 =  - V23;
                triangle(i).Vector12 = V12; 
                triangle(i).Vector21 = V21;
                triangle(i).Vector13 = V13;
                triangle(i).Vector31 = V31;
                triangle(i).Vector23 = V23;
                triangle(i).Vector32 = V32;
                triangle(i).Angle213 = (180/pi)*(atan2(norm(cross(V12,V13)),dot(V12,V13)));
                triangle(i).Angle123 =  (180/pi)*(atan2(norm(cross(V21,V23)),dot(V21,V23)));
                triangle(i).Angle132 =  (180/pi)*(atan2(norm(cross(V31,V32)),dot(V31,V32)));
                triangle(i).L12 = norm(triangle(i).Vector12);
                triangle(i).L13 = norm(triangle(i).Vector13);
                triangle(i).L23 = norm(triangle(i).Vector23);

                %Determine the boundary that each vertex is on

                   %V1
                   if ismember( v1, side1,'rows')&& ismember( v1, side2,'rows')
                         triangle(i).V1Side = 3;
                   elseif  ismember( v1, side1,'rows')
                       triangle(i).V1Side = 1;
                   elseif ismember( v1, side2,'rows')
                       triangle(i).V1Side = 2;
                   elseif ismember( v1, inner,'rows')
                       triangle(i).V1Side = 0;
                   else
                          triangle(i).V1Side = 4;
                   end

                   %V2
                   if ismember( v2, side1,'rows')&& ismember( v2, side2,'rows')
                         triangle(i).V2Side = 3;
                   elseif  ismember( v2, side1,'rows')
                       triangle(i).V2Side = 1;
                    elseif ismember( v2, side2,'rows')
                       triangle(i).V2Side = 2;
                    elseif ismember( v2, inner,'rows')
                       triangle(i).V2Side = 0;
                    else
                          triangle(i).V2Side = 4;
                    end

                   %V3
                   if ismember( v3, side1,'rows')&& ismember( v3, side2,'rows')
                         triangle(i).V3Side = 3;
                   elseif  ismember( v3, side1,'rows')
                       triangle(i).V3Side = 1;
                    elseif ismember( v3, side2,'rows')
                       triangle(i).V3Side = 2;
                    elseif ismember( v3, inner,'rows')
                       triangle(i).V3Side = 0;
                    else
                        triangle(i).V3Side = 4;
                    end

                    %Determine if the triangle is an external triangle
                     if   (triangle(i).V1Side == triangle(i).V2Side...
                           && triangle(i).V2Side == triangle(i).V3Side)...
                       ||...
                          (triangle(i).V1Side == triangle(i).V2Side...
                           && triangle(i).V3Side == 3 )...
                       ||...
                           (triangle(i).V1Side == triangle(i).V3Side...
                           && triangle(i).V2Side == 3 )...
                        ||...
                           (triangle(i).V2Side == triangle(i).V3Side...
                           && triangle(i).V1Side == 3 )...

                          triangle(i).IsExternal = 1;
                     else
                         triangle(i).IsExternal = 0;
                     end

                     %Determine if it is an end triangle
                     if triangle(i).V1Side == 3 ...
                         || triangle(i).V2Side == 3 ...
                         || triangle(i).V3Side == 3
                            triangle(i).IsEnd = 1;
                     else
                         triangle(i).IsEnd = 0;
                     end

                     %Determine the edge along the boundary for right triangles
                     if ~triangle(i).IsExternal && ~triangle(i).IsEnd
                       if triangle(i).V1Side == triangle(i).V2Side
                           triangle(i).IsRightTriangle = 1;
                           triangle(i).BoundaryEdge = [1,2];
                           if triangle(i).L13 < triangle(i).L23
                              triangle(i).RightEdge = [1,3];
                              a  = triangle(i).Vertex1;
                              b  = triangle(i).Vertex3;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                           else
                              triangle(i).RightEdge = [2,3];
                              a  = triangle(i).Vertex2;
                              b  = triangle(i).Vertex3;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                           end
                       elseif triangle(i).V1Side == triangle(i).V3Side
                           triangle(i).IsRightTriangle = 1;
                           triangle(i).BoundaryEdge = [1,3];
                           if triangle(i).L12 < triangle(i).L23
                              triangle(i).RightEdge = [1,2];
                              a  = triangle(i).Vertex1;
                              b  = triangle(i).Vertex2;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                           else
                              triangle(i).RightEdge = [3,2];
                              a  = triangle(i).Vertex3;
                              b  = triangle(i).Vertex2;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                            end
                       elseif triangle(i).V2Side == triangle(i).V3Side
                           triangle(i).BoundaryEdge = [2,3];
                           triangle(i).IsRightTriangle = 1;
                           if triangle(i).L12 < triangle(i).L13
                              triangle(i).RightEdge = [2,1];
                              a  = triangle(i).Vertex2;
                              b  = triangle(i).Vertex1;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                           else
                              triangle(i).RightEdge = [3,1];
                              a  = triangle(i).Vertex3;
                              b  = triangle(i).Vertex1;
                              triangle(i).Midpoint = [(a(1)+ b(1))/2,(a(2)+b(2))/2];
                           end 

                       else
                           triangle(i).BoundaryEdge = [0,0]; 
                           triangle(i).IsRightTriangle = 0;
                           triangle(i).Midpoint = [0,0];
                       end
                     else
                       triangle(i).BoundaryEdge = [0,0]; 
                       triangle(i).IsRightTriangle = 0;
                       triangle(i).Midpoint = [0,0];
                     end    



                if triangle(i).IsExternal == 0 && triangle(i).IsEnd == 0 

                    if    triangle(i).V1Side == 4 ...
                       || triangle(i).V2Side == 4 ...
                       || triangle(i).V3Side == 4
                            triangle(i).IsInternal = 1;
                    elseif triangle(i).V1Side ~= triangle(i).V2Side ...
                          && triangle(i).V2Side ~= triangle(i).V3Side ...
                          && triangle(i).V1Side ~= triangle(i).V3Side
                          triangle(i).IsInternal = 1;
                    else 
                        triangle(i).IsInternal = 0;
                    end
                else
                    triangle(i).IsInternal = 0;
                end

             end


% %            
% %            %Plot with external mesh removed  
% %            remove = zeros(size(tri,1),1);
% %            for i = 1:length(triangle)
% %                if triangle(i).IsExternal 
% %                    remove(i,1)  = 1;
% %                end   
% %            end
% %            index = find(remove);
% %            trimod = removerows(tri,'ind', index);
% %            figure(5)
% %            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));
% 
% 
%            %Plot with external meash and internal triangles removed
%             remove = zeros(size(tri,1),1);
%            for i = 1:length(triangle)
%                if triangle(i).IsExternal == 1 || triangle(i).IsInternal ==1
%                    remove(i,1)  = 1;
%                end   
% 
%            end
%            index = find(remove);
%            trimod = removerows(tri,'ind', index);
%            figure(6)
%            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));

            %Plot with external mesh, internal triangles, and triangles
            %sharing vertices with internal triangles.
            %Plot with external meash and internal triangles removed
            remove = zeros(size(tri,1),1);
            vertexRef = [];
           for i = 1:length(triangle)
               if triangle(i).IsExternal == 1
                   remove(i,1)  = 1;
               end
               if triangle(i).IsInternal == 1  
                   remove(i,1) = 1;
                   
                   %Remove adjacent traingels to internal ones
                   vertexRef = [vertexRef,triangle(i).Vref1];
                   vertexRef = [vertexRef,triangle(i).Vref2];
                   vertexRef = [vertexRef,triangle(i).Vref3];
                   for j = 1:size(tri,1)
                        if ismember( tri(j,1), vertexRef) || ismember( tri(j,2), vertexRef) || ismember( tri(j,3), vertexRef) 
                            remove(j,1) = 1;
                        end
                    end


               end  
               
               if triangle(i).V1Side == 4 ||triangle(i).V1Side == 4 || triangle(i).V1Side == 4
                    remove(i,1) = 1;
               end
           end

%            %Estimate the skeleton in missing sections
%             for i = 1:length(triangle)
%                if triangle(i).IsExternal == 1
%                    remove(i,1)  = 1;
%                end
%                if triangle(i).IsInternal == 1  
%                   
%                    vertexRef = [vertexRef,triangle(i).Vref1];
%                    vertexRef = [vertexRef,triangle(i).Vref2];
%                    vertexRef = [vertexRef,triangle(i).Vref3];
%                    for j = 1:size(tri,1)
%                         if ismember( tri(j,1), vertexRef) || ismember( tri(j,2), vertexRef) || ismember( tri(j,3), vertexRef) 
%                          tNormal = triangle(i).RightEdge ;
%                          tBoundary = triangle(i).BoundaryEdge;
%                          [index, loc] = ismember(tBoundary,tNormal);
%                          boundaryTerminal1 =  totalBoundary(:tNormal(loc));
%                               
%                                   tBoundary(1) a  = triangle(i).Vertex2;
%                               b  = triangle(i).Vertex1;
%                               triangle(i).Midpoint
%                                
%                         end
%                     end
% 
% 
%                end   
% 
%            end
%            index = find(remove);
%            trimod = removerows(tri,'ind', index);
%            figure(6)
%            triplot(trimod, totalBoundary(:,1), totalBoundary(:,2));


           %Plot the midpoints of the right edges
           seqnum = 0;
           midpoints = [];
           for i = 1:length(triangle)
               if triangle(i).IsRightTriangle && ~isequal(triangle(i).Midpoint,[0,0])...
                       && remove(i) ~=1
                   seqnum = seqnum + 1;
                   midpoints(seqnum,:) = triangle(i).Midpoint;
               end
           end
           figure(7)
           plot( midpoints(:,1), midpoints(:,2),'LineStyle','none','Marker','*')



%             %% Estimate of inner boundary
%            if length(B) > 1;
%                %Sample the outer boundary
%                 for i = 1:5:size(xOuter,1)
%                  xSampled = xOuter;
%                  ySampled = yOuter;
%                 end
% 
% 
%                 make_plot = 0;
%                 flag1 = 0;
%                 d = 20;
%                 [xInnerHat, yInnerHat, ~, ~, R, unv, concavity, overlap]=parallel_curve(xSampled, ySampled, d, make_plot, flag1);
% 
% 
%                 triHat = delaunay([xOuter;xInnerHat], [yOuter;yInnerHat]); 
% 
%         %         x_inner= x_inner(1:4);
%         %         y_inner = y_inner(1:4);
%         %         xx= [x_inner(1):-1:x_inner(4)];
%         %         yy = spline(x_inner,y_inner,[x_inner(1):-1:x_inner(4)]);
% 
%                 figure(8)
%                 plot( midpoints(:,1), midpoints(:,2),'LineStyle','none','Marker','*')
%                 hold on;
%                 triplot(triHat, [xOuter;xInnerHat], [yOuter;yInnerHat]);
%                 hold off
%            end

             %legend({'Curve', 'Inner Parallel', }, 'location', 'Best');

            
        
        
       
        pause;
    end


